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Abstract 

Pearl millet [Pennisetum glaucum (L.) R.Br.], a close relative of Panicoideae food crops and bioenergy grasses, offers an ideal 
system to perform functional genomics studies related to C4 photosynthesis and abiotic stress tolerance. Quantitative real- 
time reverse transcription polymerase chain reaction (qRT-PCR) provides a sensitive platform to conduct such gene 
expression analyses. However, the lack of suitable internal control reference genes for accurate transcript normalization 
during qRT-PCR analysis in pearl millet is the major limitation. Here, we conducted a comprehensive assessment of 18 
reference genes on 234 samples which included an array of different developmental tissues, hormone treatments and 
abiotic stress conditions from three genotypes to determine appropriate reference genes for accurate normalization of qRT- 
PCR data. Analyses of Ct values using Stability Index, BestKeeper, ACt, Normfinder, geNorm and RefFinder programs ranked 
PP2A, TIP41, UBC2, UBQ5 and ACT as the most reliable reference genes for accurate transcript normalization under different 
experimental conditions. Furthermore, we validated the specificity of these genes for precise quantification of relative gene 
expression and provided evidence that a combination of the best reference genes are required to obtain optimal expression 
patterns for both endogeneous genes as well as transgenes in pearl millet. 
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Introduction 

Increasing global population has raised the need of both food 
and fuel production. In addition, the growing use of fossil fuel is 
contributing to global climate changes due to elevated greenhouse 
gas emission. Pearl millet [Pennisetum glaucum (L.) R. Br., 
formerly P. americanum] is an excellent food and forage crop of 
arid to semiarid regions of the world [1,2] and a close relative of 
Panicoideae bioenergy grasses like switchgrass and foxtail millet 
[3]. It is well adapted to drought, heat, high salinity, poor soil 
fertility and low pH with an efficient C4 carbon fixation and high 
yield potential [4] . Thereby, pearl millet provides an ideal crop for 
functional genomics studies related to C4 photosynthesis and 
abiotic stress tolerance. Although several genetic engineering 
studies have been conducted in pearl millet [5,6], functional 
genomic studies under abiotic stress conditions are scanty [7]. 

Quantitative real-time polymerase chain reaction (qRT-PCR) 
provides an important platform for measuring gene expression 
changes due to its high sensitivity, specificity and wide range of 
application [8]. However, its accuracy is influenced by the 
expression stability of the internal control reference genes for 
reliable transcript normalization of target genes [9, 1 0] . An ideal 
reference gene should be constitutively and equally expressed 
across developmental stages and experimental conditions [9]. 
According to the 'golden rules' [11], identification of the most 



suitable and highly stable internal reference genes for accurate 
normalization is one of the prerequisites for qRT-PCR. So far 
most of the studies published deal with model plant species with 
known genome sequence, for e.g. Arabidopsis [12], rice [13], 
brachypodium [14]; however, relatively few studies have been 
documented in plants with limited or no genome information 
[15,16]. Thus the lack of suitable reference genes is one of the 
major limitations for gene expression studies using qRT-PCR in 
crop plants [16], including pearl millet. 

Over the past few years emphasis has been given to identify and 
validate suitable reference genes from important plant species such 
as bamboo [17], barley [18], brachypodium [14], cotton [19], 
foxtail millet [20], mustard [21], peanut [22], wheat [23,24] and 
switchgrass [25]. The commonly used traditional housekeeping 
reference genes include actin (ACT), elongation factor la, (EFla), 
glyceraldehyde-3 -phosphate dehydrogenase (GAPDH), tubulin 
(TUB), ubiquitin-conjugating enzyme (UBC) and 18S ribosomal 
RNA (18S rRNA) which are involved in basic cellular processes 
[26]. Moreover, no single traditional reference gene with stable 
constant expression across tissues and experimental conditions was 
found, thus leading to explore additional new reference genes for 
reliable normalization of qRT-PCR data [26]. Recent reports 
illustrated that F-box/kelch-repeat protein (F-box), phosphoenolpyr- 
uvate carboxylase-related kinase (PEPKR), protein phosphatase 2A 
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(PP2A) and TIP41 -like family protein (TIP41) genes are superior 
compared to traditional reference genes [17,21,22,27]. Several 
statistical algorithms, namely, Stability Index [28], ACt [29], 
BestKeeper [30], geNorm [31], NormFinder [32], and RefFinder 
[33] have been employed for proper validation and stability 
ranking of the best reference genes for qRT-PCR data normal- 
ization in numerous plant species. However, to the best of our 
knowledge, no systematic analysis for the selection of suitable 
reference genes for qRT-PCR analysis in pearl millet has been 
reported. Therefore, a comprehensive validation of reference 
genes under different experimental conditions for accurate 
transcript normalization is needed in pearl millet. 

In this work, we evaluated 18 potential candidate reference 
genes in 234 samples from three important pearl millet genotypes 
using qRT-PCR. Expression patterns of these genes were 
monitored in tissue samples under different developmental 
processes, hormone treatments and abiotic stress conditions. 
Expression stability of these genes was validated using six statistical 
algorithms in order to assign appropriate reference genes suitable 
to each experimental condition for accurate transcript normaliza- 
tion. Our results showed that sets of genes are appropriate for 
accurate transcript quantification of endogenous genes as well as 
transgenes from different tissue samples. We further illustrated 
detailed expression patterns of three essential pearl millet 
endogenous genes specific to development, hormonal stimuli and 
abiotic stresses. 

Materials and Methods 

Plant materials 

Pearl millet (Pennisetum glaucum [L.] R. Br.) genotypes 
ICMR01004, IPCI1466 and IP300088 were used in this study. 
Seeds of ICMR01004 and IPCI1466 were obtained from the 
International Crops Research Institute for the Semi-Arid-Tropics 
(ICRISAT), India, while seeds of IP300088 were acquired from 
the Germplasm Resources Information Network's (GRIN), USA. 
Seeds were kept in wide mouth polypropylene botdes (VWR) and 
stored in a seed vault at 9°C with a relative humidity of 50%. 

Developmental tissue samples 

For developmental tissue samples, three genotypes were grown 
in 5 liter pots containing agronomy mix (equal parts of redwood 
compost, sand and peat moss) under greenhouse conditions of 
16 h day/8 h night photoperiod at 30±2°C until maturity. Plants 
were watered every alternate day with tap water and fertilized 
biweekly. Tissue samples of vegetative and reproductive stages 
included callus 30DPC (days post culture), leaf 7DPS (days post 
sowing), leaf 15DPS, leaf 30DPS, node, internode, sheath, flag 
leaf, panicle, peduncle and root of 60DPS plant, and 30DPH (days 
post harvest) seeds. A total of 108 tissues samples comprising of 12 
vegetative and reproductive stages from three genotypes in three 
biological replicates were harvested by immediate quick freezing in 
liquid nitrogen in 2 ml SealRite microcentrifuge tubes. 

Hormone treatments 

Seeds of 30DPH were soaked in 70% (v/v) ethanol for 1 min 
followed by washing in 2.5% (v/v) sodium hypochlorite solution 
containing 0.1% (v/ v) Tween 20 for 1 5 min and rinsed thoroughly 
with sterile distilled water. Surface sterilized seeds were grown in 
PhytoCon culture vessels (Phytotechnology Laboratories, Over- 
land Park, KS, USA) containing half strength Murashige and 
Skoog (MS) medium for 14 days. Seedlings were kept in sucrose 
free liquid half strength MS medium for 24 h. Seedlings of 15DPG 
(days post germination) were transferred to PhytoCon culture 



vessels (Phytotechnology Laboratories) containing liquid half 
strength MS supplemented with 100 uM abscisic acid (ABA, 
Sigma, St. Louis, MO, USA), 50 uM brassinolide (Bra, Sigma), 
50 |jM gibberellic acid (GA, Sigma), 50 uM indole-3-acetic acid 
(IAA, Sigma), 100 uM methyl jasmonate (Meja, Sigma), 100 uM 
salicylic acid (SA, Sigma), 100 U.M Zeatin (Zea, Sigma) and 
incubated for 6 h. Leaves from a total of 72 samples from seven 
treatments in three biological replicates including one untreated 
control of three genotypes were harvested and immediately frozen 
as mentioned in the earlier section. 

Abiotic stress conditions 

In the dehydration stress treatments, seedlings of 1 5DPG (same 
as hormone treatments) were kept in 400 |xM mannitol solution 
for 6 h. For drought and salinity stresses, water supply was 
withheld and 300 mM sodium chloride (NaCl) solution was 
provided for 5 days to 30DPS plants, respectively. Cold and heat 
stresses were carried out by maintaining 30DPS plants at 4± 1°C 
and 42±1°C, respectively for 6 h for 3 consecutive days. Stress 
symptoms were monitored visually by the appearance of leaf 
rolling and yellowing, as well as by measuring stomatal conduc- 
tance and photosynthesis rates of plants using a LI-COR 6400-40 
with an integrated fluorescence chamber head (LI-COR, Lincoln, 
Nebraska, USA) after the stress treatments. 

Candidate reference genes selection and primer design 

Locus identifiers (IDs) of Arabidopsis and rice potential 
candidate reference genes were obtained from previously pub- 
lished work (Table 1). Orthologous locus IDs from foxtail millet 
(Setaria italica) were identified using locus search from Phyto- 
zome. GenBank accession numbers were obtained from National 
Center for Biotechnology Information (NCBI) using BLASTN. 

A total of eighteen genes were chosen for primer design using 
Primer3Plus software (http://www.bioinformatics.nl/cgi-bin/ 
primer 3plus/primer3plus.cgi) [34] considering the parameters 
specific for qRT-PCR. The sequences with detailed parameters 
for each primer pair are given in Table SI. 

RNA isolation and cDNA synthesis 

A total of 100 mg of frozen plant material was ground to fine 
powder in a 2 ml SealRite microcentrifuge tube using 3.2 mm 
stainless steel beads and an automated shaker SO-10M (Fluid 
Management, Wheeling, IL, USA). Total RNA was isolated from 
plant samples using the RNeasy plant mini kit (Qiagen, Valencia, 
CA, USA) according to the manufacturer's procedure. A first set of 
on-column DNAse I (Qiagen) digestion was carried out during the 
RNA extraction steps. The integrity of RNA samples were 
checked by 1% (w/v) agarose gel. The quantity and quality of 
RNA samples were also checked using a NanoDrop ND-1000 
(NanoDrop Technologies, Wilmington, DE, USA). RNA samples 
with 260/280 ratio between 1.9 to 2.2 and 260/230 ratio between 
2.0-2.5 were used for cDNA synthesis. To completely eliminate 
DNA contamination, 1 |xg of total RNA was subjected to gDNA 
wipeout reaction using the QuantiTect reverse transcription kit 
(Qiagen) followed by first strand cDNA synthesis in a 20 ul 
reaction mixture using an optimized blend of oligo-dT and 
random primers according to manufacturer's instructions and 
stored at -20°C. 

PCR and qRT-PCR 

Specific amplification from cDNA was checked by PCR 
following the protocol described earlier [35] using 1 n.1 of cDNA, 
1 0 mM clNTPs, 1 U.M each of forward and reverse primers and 
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Table 1. Expression level of the selected candidate reference genes tested in pearl millet. 



Genes 


Description 


Arabidopsis 


Rice 




Foxtail millet/Pearl millet 3 


Ct b ±SD 


CV±SD 


ACT 




Ml I gzzozu 


LOC 


_usuDg jozyu 


jiuzz.3/ zm.g/niviz i f jduu 


28 2~ l ~2 7 


1 8 + o 5 


rvr 


Cyclophilin, peptidyl-prolyl cis-trans isomerase 


AlDg 3D I uu 


i r\r 


_usut>g i yo i u 


diu 1 om.g 


31 .5 it 3.0 


3.3 ± 1 .0 


eEFlot 


Eukaryotic elongation factor 1 alpha 


AiDgou^yu 


LUL_ 


.usuoguoiou 


jiuzzujym.g/troy'4- 1 dd 


24.1 it 3. 2 


1 .8 ±0.5 


PQY 
rtSA 


F-box domain containing protein 


AlDg I Dl I u 


1 C\C 


_usu^gD / zyu 


jiuzz i jom.g 


25.3 it2. 8 


2.6 ±0.7 


CAonu 
KjArUrt 


Glyceraldehyde-3-phosphate dehydrogenase 


Aijgu^ i zu 


i r\r 


.usuoguozyu 


diu i^u^^m.g/utjjycs i u/ 


22.8it3.1 


2.6±0.6 


elF4o2 


Eukaryotic initiation factor 4a2 


At-1 rtKA ~)~7f\ 
AI I gD4Z/U 


M 
IN 




JlUUDD4om.g/ tUODDD 3D 


23.2it2.1 


2.9±0.7 


rtrl\ri 


Phosphoenolpyruvate carboxylase kinase related 


AI I g I ZDoU 


1 C\C 


.usuoguoooz 


DlUUdz/ jm.g/rno/z/oo 


25.6 it 1 .5 


1 .3 — 0.4 


PP2A 


Protein phosphatase 2A 


AI I g I LKf 3\J 


LOC 


_usuzg i z dou 


<;ini 7QQ1m n 

jiu i / oyzm.g 


25 6~ l ~2 5 


12 + 03 


RCA 


Rubisco activase 


AlZg 3¥ 1 3\J 


LOC 


_us i i g^f/y/u 


jiuzo^f 1 4m.g 


24 9~ l ~2 9 


2 4+0 6 


DMIVIUC 


S-adenosyl methionine decarboxylase 




LOC 


.usu^+g^+zuyu 


<;ini rnftim n 
jiu i uztszm.g 


26 2~ l ~5 2 


4 8 + 1 2 


TUA 


Tubulin alpha 


At1g04820 


LOC. 


_Os03g51600 


Si035654m.g 


23.1 ±3.4 


2.1 ±0.5 


TIP41 


Tonoplast intrinsic protein 


At4g34270 


LOC. 


_Os03g55270 


Si036884m.g 


28.5 ±1.5 


1.1 ±0.3 


UBC2 


Ubiquitin-conjugating enzyme 2 


At5g25760 


LOC. 


_Os02g42314 


Si018564m.g 


29.8±3.1 


1.7±0.5 


UBC18 


Ubiquitin-conjugating enzyme 18 


At5g42990 


LOC. 


_Os12g44000 


Si023498m.g 


26.3±2.2 


2.7 ±0.6 


UBQ5 


Ubiquitin 5 


At2g47110 


LOC. 


_Os01g 22490 


Si003209m.g 


23.5±2.1 


1.3 ±0.4 


UNK 


Transmembrane protein 56 


At1g31300 


LOC. 


_Os01g56230 


Si002525m.g 


27.9±1.7 


2.6 ±0.6 


18S rRNA 


18S ribosomal RNA 


N 


N 




KC201690 


24.0±4.9 


5.7±1.3 


25S rRNA 


25S ribosomal RNA 


N 


N 




AB197128 


9.1 ±1.8 


3.6±0.3 



a Locus identifiers of selected candidate reference genes for foxtail millet and/or GenBank accession numbers for pearl millet with orthologous from Arabidopsis and 
rice are listed. 

b The expression levels of the candidate genes obtained during qRT-PCR experiments of total samples (n = 234) are presented as mean threshold cycle (Ct) values. SD, 
standard deviation; CV, coefficient of variance; N, no corresponding locus identifier or accession number. 
doi:1 0.1 371 /joumal.pone.01 06308.t001 



one unit Taq polymerase in a 10 (J.1 total reaction mixture. The 
amplification program was as follows: 5 min at 95°C; followed by 
30 cycles of 95°C for 30 sec, 58°C for 15 sec, 72°C for 30 sec; and 
a final extension of 72°C for 10 min followed by electrophoresis on 
3% (w/v) agarose gel. 

For qRT-PCR, cDNAs were diluted to 20 times into a final 
volume of 400 |J.l and the reactions were performed as described 
previously [36] in an optical 96 well plate (Applied Biosystems, 
Foster City, CA, USA) containing 1 |xl of diluted cDNA, 200 nM 
of each gene specific primer and 2.5 |il of 2X Fast SYBR Green 
PGR master mix in a 5 ul total volume using a StepOnePlus™ 
real time PCR system (Applied Biosystems) equipment. The qRT- 
PCR reactions were conducted following the fast thermal cycles: 
50°C for 2 min, 95°C for 20 sec, followed by 40 cycles of 95°C for 
3 sec and 60°C for 30 sec. After 40 cycles, the specificity of the 
amplifications was tested by heating from 60°C to 95°C with a 
ramp speed of 1.9°C/min, resulting in melting curves. The 
threshold cycle (Ct) value was automatically determined for each 
reaction by the real time PCR system with default parameters. 
Raw data (not baseline corrected) of fluorescence levels and the 
specificity of the amplicons were checked by qRT-PCR dissoci- 
ation curve analysis using StepOne Software (v2.3). The baseline 
correction and linear regression analysis on each amplification 
curve including the efficiencies (E) of the polymerase chain 
reactions were calculated based on the slope of the line 
(E— 10 sope ), considering an ideal value range (1.8^ £&2) and 
correlation (R 2 >0.9) using the LinRegPCR software [37] The 
final Ct values were the mean of three biological replicates and the 
coefficient of variance (CV) was calculated to evaluate the 
variation of Ct values for each gene. Each qRT-PCR reaction 
set included water as a negative no-template control (NTC) instead 
of cDNA. 



Analysis for expression stability of reference genes 

Five different types of computer-based programs, Stability 
Index [28], delta (A)Ct [29], BestKeeper [30], geNorm [31] and 
NormFinder [32] methods were used to rank and compare the 
stability of candidate reference genes across all the experimental 
sets. For Stability Index, ACt, BestKeeper programs, the Ct value 
for each candidate reference gene was used to determine its 
relative expression stability. For NormFinder and geNorm, relative 
expression values were calculated from 2 AACt using the formula 
applied before [31]. Overall recommended comprehensive 
geomean ranking values of the best reference genes were obtained 
using the ranking results of four algorithms, except Stability Index, 
in RefFinder [33]. The pairwise variation (Vn/Vn+1) between two 
sequential normalization factors (NFn and NFn+1) were estimated 
using geNorm software provided in qBasePlus (v2.4) [38] package 
for best and minimal number of reference genes needed to 
calculate an optimal normalization. 

Validation of reference genes 

Six genes were chosen to determine their differential expression 
after accurate normalization across five experimental sets using 
single and/or best combinations of reference genes (Table S2). 
Primer design and qRT-PCR reactions were followed as 
mentioned before. The average Ct value was calculated from 
three biological replicates and used for relative expression 
analyses. Normalization of the gene of interest in developmental 
tissue samples was calculated using the ACt values as previously 
described [12], while relative expression of genes of interest in 
hormone treatments and abiotic stress conditions was measured as 
suggested before [39]. The expression fold change value was 
represented as relative expression (2 AACt ). Statistical significant 
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differences in gene expression patterns were evaluated using 
Tukey's range test in JMP (v7.0.2). 

Transformation of pearl millet 

Particle bombardment-mediated transformation of immature 
zygotic embryo derived calli was carried out using PDS-1000 He 
biolistic device (Bio-Rad, Hercules, CA) following the protocol 
described earlier [6]. Zygotic embryos were isolated from surface 
sterilized seeds and cultured on MS medium supplemented with 
2,4-D (2.5 mg/1), maltose (30 g/1), pH 5.8 for callus formation. 
Particle bombardments were conducted using pCAMBIA1201 
and pCAMBIA1302 vectors plasmid DNA (250 ng/shot) precip- 
itated onto 0.6 u.m gold particles (Bio-Rad) at a helium pressure of 
1,100 psi following the protocol describe previously [6]. Expres- 
sion of ^-glucuronidase (gus) reporter gene was performed as 
mentioned earlier [40], while green fluorescent protein (gfp) 
expression was monitored using a fluorescence stereomicroscope 
(Leica MZ FLIII) coupled with a SPOT Insight CCD camera. 

Results 

Identification of candidate reference genes 

We found locus identifiers and/ or GenBank accession numbers 
of selected Arabidopsis and rice candidate reference genes from 
previous published work (Table 1). We identified orthologous 
locus IDs and/ or GenBank accession numbers of these potential 
candidate reference genes from foxtail millet, a close relative of 
pearl millet, using orthologous group search in Phytozome and/ or 
BLASTN search in NCBI GenBank. We selected a total of 18 
genes for accurate transcripts normalization during gene expres- 
sion study using qRT-PCR in pearl millet. These genes included 
both traditional housekeeping as well as several new reference 
genes namely, actin (ACT), cyclophilin (CYC), eukaryotic elonga- 
tion factor 1 alpha (eEFla), F-box domain containing protein 
(FBX), glyceraldehyde-3 -phosphate dehydrogenase (GAPDH), eu- 
karyotic initiation factor 4a2 (eIF4a2), phosphoenolpyruvate 
carboxylase-related kinase (PEPKR), protein phosphatase 2 A 
(PP2A), rubisco activase (RCA), S-adenosyl methionine decarbox- 
ylase (SAMDc), alpha tubulin (TUA), tonoplast intrinsic protein 
(TIP41), ubiquitin-conjugating enzyme 2 (UBC2), ubiquitin- 
conjugating enzyme 18 (UBC18), ubiquitin 5 (UBQ5), transmem- 
brane protein 56 [LINK), 18S ribosomal RNA (18S rRNA) and 
25S ribosomal RNA (258 rRNA) (Table 1). 

Due to insufficient availability of sequence information in the 
NCBI GenBank, in addition to the sequences of pearl millet 
obtained from Genbank, we used full length transcript sequences 
from the foxtail millet to design the gene specific primers for qRT- 
PCR. Primer pairs were designed to anneal near the 3' end or at 
the 3' UTR of each gene using the Primer3Plus software following 
the parameters: length: 20 ±3 mer; product size range: 50-200 
base pair; melting temperature: 60°C±3°, guanine-cytosine (GC) 
content: ~50% including absence for hairpin structures, self- 
dimers and weak or no self-complementarities at the 3' end (Table 
SI). 

Sample size, RNA quality and qRT-PCR conditions 

We tested the expression of these potential candidate reference 
genes and quantified the Ct values using qRT-PCR in a total 
experimental set of 234 samples (Table 1). These included 
developmental tissues, hormone treatments and abiotic stress 
conditions from three pearl millet genotypes ICMR01004, 
IPCI1466 and IP300088 (Table 2). The developmental tissues 
experimental set included 108 samples from 12 vegetative and 
reproductive stages [callus 30DPC, seed 30DPH, leaf 7DPS, leaf 



15DPS, leaf 30DPS, and node, internode, sheath, flag leaf, 
panicle, peduncle and root from 60DPS plants] , whereas hormone 
treatments and abiotic stress conditions included 72 and 54 
samples from 8 [control (without treatment), ABA, Bra, GA, IAA, 
Meja, SA and Zea] and 6 [control (without stress), dehydration 
(mannitol), drought (no water), cold, heat and salinity] sets of 
samples (Tables 2, S3 to S5), respectively. The fifth experimental 
set comprised of Ct values from 78 tissue samples from each of the 
three pearl millet genotypes (Table 2). We isolated high quantity 
[368.7±63.3 ng/u.1 (mean ± standard deviation, SD where 
n = 234)] and quality (average 260/280 ratio of 2.0±0.1 and 
260/230 ratio of 2.1 ±1.6) of total RNA using the guanidinium 
thiocyanate-based RNeasy plant mini kit. The complete absence 
of DNA contamination was confirmed by qRT-PCR after two 
steps of DNAse treatments (first on-column and second gDNA 
wipeout reaction) for each sample. The reverse transcriptase 
reactions were primed using an optimized blend of oligo-dT and 
random primers provided in the kit in order to amplify transcripts 
from both highly and weakly expressed genes. 

Accuracy and efficiency of amplification 

To determine the accuracy of primers designed in this study to 
specifically amplify potential target candidate reference genes we 
performed PCR and qRT-PCR using either crude and/ or diluted 
cDNAs. We obtained a single amplified product of the expected 
size in agarose gel electrophoresis (Figure SI) and the presence of 
one dominant peak of the specific amplicon in melt curve analysis 
(Figure S2), respectively. Further, a two-step qRT-PCR protocol 
for cDNA synthesis and cDNA amplification in successive steps 
reduced the undesired primer dimer formation using SYBR 
Green. No detectable amplifications in the no-template controls 
(NTCs) confirmed the absence of primer dimers or non-specific 
products (Figure S2). We determined the PCR efficiency (E) of 
each primer pair from the amplification plots of all amplification 
profiles using LinRegPCR software. The mean E values with SD 
for all primer pairs across all experimental samples from three 
biological replicates are given in Table S 1 . Primer pairs of most of 
the genes exhibited no significant differences in E values and 
displayed PCR efficiencies of more than 1.90, while primer pair 
for CYC and 258 rRNA showed PCR efficiencies of 1.87±0.03 
and 1.85±0.04 (Table SI). We further calculated correlation 
coefficients (R 2 ) of PCR efficiency values to evaluate the 
amplification curves. Except TAU (R 2 = 0.89), the rest of the 
primer pairs revealed i? 2 >0.90 from all reactions (Table SI). 

Expression levels of candidate reference genes 

Expression levels of all the candidate reference genes were 
measured by monitoring the Ct values in the qRT-PCR reactions. 
We analyzed all the Ct values under five groups which included 
total [first experimental set (n = 234), Table 1], developmental 
tissues [second experimental set (n= 108), Tables 2 and S3], 
hormone treatments [third experimental set (n= 72), Tables 2 and 
S4], abiotic stress conditions [fourth experimental set (n = 54), 
Tables 2 and S5] and genotypes [fifth experimental set (n = 78), 
Table 2] . In the first total experimental set the mean Ct values of 
the 18 candidate reference genes revealed a minimum of 9.1 ±1.8 
and a maximum of 3 1.5 ±3.0 for highest and lowest expression 
levels for 25S rRNA and CYC genes, respectively, while most of 
the values were distributed between 22.8±3.1 to 29.8±3.1 
(Table 1). The mean Ct values of SAMDc (26.2±5.2) with highest 
SD indicated less stability as compared to TIP41 (28. 5 ±1.5) with 
lowest SD showing relatively stable expression in the total 
experimental set (Table 1). Similarly, we noticed large SDs of Ct 
values for SAMDc and 18S rRNA indicating a more variable 
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Table 2. Expression levels of candidate reference genes across four experimental sets of pearl millet. 



Genes 


Developmental tissues Ct±SD* 


Hormone treatments Ct±SD 


Abiotic stresses Ct±SD 


Genotypes Ct±SD 


ACT 


28 1 _h 2 8 


27 6 + 1 7 


29 3 3 5 


28 3 "*"0 7 


CYC 


29 8 + 3 2 


33 2 -1- 1 1 


32 3 ~*~3 0 


31 7 + 1 4 


eEFIoL 


23 0^3 6 


25 3 -1- 1 9 


24 3 —i— 3 1 


24 2~^0 9 


FBX 


25 1 + 3 7 


24 5 -1- 1 3 


26 8~ , ~2 2 


25 4~ l ~0 9 


GAPDH 


23 1 ~ 1 ~3 4 


21 6~*" 1 2 


24 O - ^ 6 


22 9"*" 1 0 


elF4a2 


23 0 + 2 1 


22 5 -1- 1 1 


24 7~ l ~2 4 


23 3 1 0 


PEPKR 


25 3 + 1 3 


25 4~ l ~0 7 


26 3 ~ t_ 2 3 


25 6 + 0 5 


PP2A 


25 3 + 2 9 


24 8 + 1 2 


26 5 + 2 6 


25 7 + 0 7 


RCA 


25 8 + 3 5 


24 1 + 0 9 


24 5 + 3 0 


24 8 + 0 7 


SAMDc 


27 o + 5 5 


24 9+4 0 


26 6 + 5 8 


26 2 + 0 9 


TUA 


22.8±3.9 


21.4±1.0 


26.0±2.8 


23.4±1.9 


TIP41 


28.3 ±1.7 


28.1 ±0.7 


29.3 ±1.7 


28.6 ±0.5 


UBC2 


29.9±3.3 


28.8±2.0 


30.8±3.5 


29.8±0.8 


UBC18 


26.3±2.7 


25.8±0.7 


26.9±2.4 


26.3 ±0.4 


UBQ5 


23.7±2.5 


22.3 ±0.8 


24.7±2.0 


23.5 ±1.0 


UNK 


27.7±2.0 


27.6±1.1 


28.7±1.7 


28.0 ±0.5 


18S rRNA 


22.5±5.4 


26.4±3.0 


23.4±5.0 


24.1 ±1.7 


25S rRNA 


9.0±2.3 


8.7±0.6 


9.8 ±1.9 


9.1 ±0.5 



*, Data are represented as mean threshold cycle (Ct) values from all analyzed samples in each individual experimental set with standard deviation (SD). 
doi:1 0.1 371 /journal.pone.01 06308.t002 



expression in the other four experimental sets, while little variation 
of Ct values was detected for rest of the genes (Table 2). In the 
second experimental set, high Ct values suggested low expression 
of all the candidate reference genes in the seeds compared to other 
developmental tissues (Table S3). The third and fourth experi- 
mental sets revealed elevated Ct values of these genes under SA 
treatment (Table S4), and heat and salinity stress conditions (Table 
S5) as compared to their respective controls. Furthermore, we 
evaluated the expression levels of candidate reference genes by 
calculating the CV of the Ct values. Among the four experimental 
sets, TIP41 showed the lowest CV value (1.1 ±0.3), while SAMDc 
and 18S rRNA revealed a greater variation in expression levels 
due to their high CV values (4.8±1.2 and 5.7± 1.3, Table 1). 

Stability ranking of the candidate reference genes 

We used Stability Index (SI), BestKeeper, ACt, NormFinder, 
geNorm and RefFinder programs to identify the best reference 
genes for qRT-PCR data normalization in pearl millet. These 
programs allowed us to establish a stability ranking of each 
candidate reference gene using Ct values across the experimental 
sets and condition-specific levels (Tables 3-7). 

The SI was calculated from the multiplication of the slope 
(a value of the regression analysis of geomeans and overall means) 
with CV considering the fact that gene with lowest SI from low 
slope and low CV provided the best reference gene. In the first 
total experimental set, PEPKR, PP2A and TIP41 with SI values 
of 0.05, 0.06 and 0.09 were the top three candidates, respectively, 
whereas 18S rRNA with highest SI value of 1.97 was the least 
preferred choice of reference gene (Table 3). Based on SI values, 
PEPKR (SI of 0.06) and TIP41 (SI of 0.16) were the two best 
candidates, while SAMDc (SI of 2. 40) was the worst candidate for 
normalization of gene expression in developmental tissue samples 
(Table 4). Analysis of SI values of reference genes in the third and 
fourth experimental sets revealed TIP41 as the superior candidate 



with the smallest SI of 0.02 and 0.18, respectively for transcript 
normalization under hormone treatments and abiotic stress 
conditions (Tables 5 and 6). Among the three genotypes of pearl 
millet, TIP41 (SI of 0.16) was the top ranked reference gene for 
normalization of gene expression (Table 7). 

The BestKeeper program determines the stability ranking of the 
reference genes based on the percentage of crossing point (%CP) 
to the BestKeeper Index and the SD from the geometric mean of 
the candidate reference genes Ct values, where the genes with 
lowest CP and SD values are identified as the best reference genes 
for normalization. In this study, BestKeeper analyses of the total 
experimental samples identified PEPKR (2.72±0.69), TIP41 
(3.16±0.90) and PP2A (3.24±0.90) with lowest CP±SD values 
(Table 3), where genes with SD<1 are considered as stable. In the 
developmental tissues, hormone treated, abiotic stressed and 
genotype experimental sets many genes showed SD<1, while 
the most stable reference genes were PEPKR (2.78±0.70), TUA 
(0.70±0.23), TIP41 (1.60±0.47) and ACT (1.09±0.29), respec- 
tively (Tables 4-7). 

The ACt method compared the relative expression of a 
reference gene with other candidate reference genes within each 
sample, thereby ranked the genes based on the average of STDEV 
or SD. Analyses using this program exhibited PP2A, UBC2, 
TIP41, UBQ5 and TIP41 with average STDEV values of 1.32, 
1.33, 0.56, 0.97 and 0.64 as the most suitable reference genes for 
normalization in total, developmental, hormone treated, abiotic 
stress and genotypes experimental sets of pearl millet, respectively 
(Tables 3-7). 

NormFinder ranks all candidate reference genes based on intra- 
and inter-group variations of expression stabilities by measuring 
the stability value (SV) for each reference gene. In our study, the 
NormFinder identified PP2A, TIP 41 and PEPKR with SV of 
0.43, 0.57 and 0.61, as the top three optimal reference genes for 
transcript normalization in the total tissue samples (Table 3). The 
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0.30 




■ Total 

hi Developmental tissues 
m Hormone treatments 
□ Abiotic stresses 
b Genotypes 



Figure 1. Estimation of pairwise variation to determine the optimal number of control reference genes required for accurate 
normalization using geNorm. Pairwise variation (V, Vn/Vn+1) was calculated between successively ranked normalization factors NFn and NFn+1. 
Arrowheads on the bar graph indicate the minimum number of genes required at the cut-off value 0.15 [31]. The V between the normalization factors 
of the two first-ranked and the three first-ranked is represented by V2/3 and so on, respectively. 
doi:1 0.1 371 /journal.pone.01 06308.g001 



NormFinder analyses in the developmental tissues (SV of 0.31), 
hormone treatments (SV of 0.13), abiotic stress conditions (SV of 
0.28) and genotypes (SV of 0.03) experimental sets of pearl millet 
recognized TIP41 as the most suitable reference gene (Tables 4-7). 

In addition, we also examined the stability ranking of candidate 
reference genes using geNorm program (Tables 3-7). The 
geNorm statistical algorithm determines the normalization value 
(MV) based on the geometric mean of multiple reference genes 
and mean pair-wise variation of a gene from all other reference 
genes in each set of samples. In both first and second experimental 
sets, the two best reference genes were PP2A \ TIP41 with the 
lowest MV of 0.46 and 0.32, whereas UBC2 with MV of 0.49 and 
0.36 remained the third most suitable gene for transcript 
normalization in total and developmental tissues, respectively, as 
determined by the geNorm (Tables 3-4). The most preferred 
genes for normalization in hormone treatments and abiotic stress 
conditions were TIP41 \ UBQ5 (MV of 0.16) and PP2A \ TIP41 
(MV of 0.39), respectively (Tables 5-6), while TIP41 \ACT had 
the lowest MV of 0.05 in the genotypes of pearl millet (Table 7). In 
addition, geNorm analyses revealed significantly high stability of 
several reference genes with MV of less than the cut-off range of 
1.5 (Tables 3-7). 

We further compared all the data generated by SI, BestKeeper, 
ACt, NormFinder and geNorm programs using recommended 
comprehensive ranking method in RefFinder software to confirm 
the stability ranking of reference genes for accurate transcript 
normalization across the experimental sets (Tables 3-7). The 
overall ranking of the best reference genes in total and categorized 
experimental sets according to RefFinder are given in Tables 3-8. 

We next applied the geNorm software to calculate the Vn/Vn+1 
between NFn and NFn+1 to determine the best combination of 
reference genes required for precise transcript quantification across 
different sets of experiments. Figure 1 summarizes the V values 
from the combination of reference genes and shows that a number 
of genes are required for reliable normalization of gene expression 
data among different experimental sets (Table 8). 

Accurate normalization of gene and transgene 
expression using optimal combination of reference 
genes 

In order to validate the selection of the best reference genes for 
accurate normalization of gene expression, we chose PEPC 
(phosphoenolpyruvate carboxylase), ERF (ethylene response 



factor) and DREB (dehydration responsive element binding) genes 
to determine the relative transcript levels using qRT-PCR (Table 
S2). We monitored the expression of PEPC, an essential gene for 
C4 photosynthesis, in developmental tissue samples, whereas the 
expression pattern of two transcription factors, ERF and DREB, 
known to be regulated during abiotic and biotic stresses, were 
examined in hormone treated and abiotic stressed samples. 
Relative transcript levels of these genes were calculated after 
normalizing with the best ranked candidate reference genes 
as determined by geNorm and recommended by RefFinder 
(Table 8). Transcript abundance of PEPC when normalized using 
single top ranked reference genes, PP2A, TIP41 and UBC2, 
revealed bias effect on the relative expression patterns (Figure 2). 
Furthermore, transcript normalization using a combination of two 
{PP2A+TIP41) and three (PP2A+TIP41+UBC2) reference genes 
showed much stable and constant expression profiles across tissues 
(Figure 2). Similarly, relative expression patterns of ERF and 
DREB in hormone treatments and abiotic stress conditions were 
affected by the selection of the reference gene or combination of 
genes, respectively (Figures 3-4). As predicted, a strong bias in the 
relative expression pattern of PEPC, ERF and DREB was 
obtained when the least stable gene was used for normalization. 
Incorporation of TIP41 and UBC2 or UBQ5 during expression 
analyses neutralized the unwanted changes of transcript abun- 
dance to allow accurate normalization of PEPC, ERF and DREB. 
Overall expression of PEPC was significandy high in flag leaf and 
sheath as compared to nodal tissues of pearl millet genotypes 
(Figure 2). In the hormone treatments experimental set, Zea 
enhanced 2-fold expression of ERF in pearl millet genotype 
ICMT01004 and IPCI1466 compared to other hormones tested 
(Figure 3). The expression of DREB was up-regulated during 
drought followed by heat stresses in all the three genotypes 
(Figure 4). Genotypes showed differential expression patterns of 
these genes as well (Figures 2-4). 

We also monitored the transcript abundance pattern of 
^-glucuronidase (gus), green fluorescent protein (gfp) and hygro- 
mycin phosphotransferase (hpt) expressing transgenes in transgenic 
pearl millet calli. Calli of three pearl millet genotypes were 
bombarded with CaMV35S::gus (pCAMBIA1201) and 
CaMV35S::gfp (pCAMBIA1302) constructs and transient expres- 
sion of both gus and gfp reporter genes were visualized after 5 days 
(Figure S3). Expressions of gus, gfp and hpt genes were examined 
in transformed calli selected on hygromycin (30 mg/1) after 30 
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Table 8. Summary of the best combination of reference genes for accurate normalization across five experimental sets of pearl 
millet using geNorm and RefFinder programs. 





Experimental sets 


Total 


Development tissues 


Hormone treatments 


Abiotic stresses 


Genotypes 


Best combination 


V2/3 


V3/4 


V2/3 


V3/4 


V2/3 


Pairwise variation (V) a 


0.070 


0.138 


0.056 


0.090 


0.110 


Reference control genes 


PP2A 


PP2A 


TIP41 


PP2A 


TIP41 




TIP41 


TIP41 


UBQ5 


TIP41 


ACT 






UBC2 




UBQ5 





a Pairwise variation (V) represents the optimal combination of reference control genes required to pass the suggested cut-off value 0.15 [31]. A single common reference 
control gene for expression study across experimental sets is highlighted in gray. 
doi:1 0.1 371 /journal.pone.01 06308.t008 
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Figure 2. Validation of PEPC gene expression after normaliza- 
tion using optimal number of control reference genes in 
developmental tissue samples from genotypes (A) ICMR01004, 
(B) IPCI1466 and (C) IP300088. Results are presented as mean 
relative expression with SD from three biological replicates after 
normalization using the best combination of reference genes recom- 
mended by geNorm and RefFinder (see Table 8) for developmental 
tissue samples. Leaf- 7D, 15D and 30D represent 7DPS, 15DPS and 
30DPS leaf samples while flag leaf, sheath, node and internode are from 
60DPS plants. Different letters on the bars indicate significant 
differences at the P<0.05 level as tested by Tukey's Range Test. 
doi:1 0.1 371 /journal.pone.01 06308.g002 
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Figure 3. Validation of ERF gene expression after normalization 
using optimal number of control reference genes in hormone 
treated samples from genotypes (A) ICMR01004, (B) IPCI1466 
and (C) IP300088. Data are presented as mean relative expression 
with SD from three biological replicates after normalization using the 
best combination of reference genes recommended by geNorm and 
RefFinder (see Table 8) for hormone treatments. ABA (abscisic acid), Bra 
(brassinolide), GA (gibberellic acid), IAA (indole-3-acetic acid), MeJa 
(methyl jasmonate), SA (salicylic acid) and Zea (zeatin) treatments of 
15DPG plants. Different letters on the bars indicate significant 
differences at the PsO.05 level as tested by Tukey's Range Test. 
doi:10.1 371/journal.pone.01 06308.g003 



PLOS ONE | www.plosone.org 



12 



August 2014 | Volume 9 | Issue 8 | e106308 



Assessing Reference Genes for qRT-PCR in Pearl Millet 



^ 60 0 
c 50 0 



> 20 0 
J5 

£ 10 0 

0.0 

g 1600 0 
1400.0 

c 

.2 1200.0 
10 

<2 1000 0 

X 800 0 

(1) 600 0 
> 

ro 400 o 

Ct 200 0 

0.0 

Q 300 0 

250.0 

c 
O 

'« 200 0 

£ 

Q- 1500 
0 



| 500 
DC 



1 



DPP2A 

OTIP41 

■ U805 

■PP2A+TP4) 

mPP2A*TIP41*UBQ5 



b b b 




a a a a a 



[*] b b b b 



d d 

no 



b 

li 



Hi 



T a a a a 

iViiii 



Onl A 



a a a a a 



Dehydration Drought 



Heat 



Cold 



Figure 4. Validation of DREB gene expression after normaliza- 
tion using optimal number of control reference genes in 
genotypes (A) ICMR01004, (B) IPCI1466 and (C) IP300088 
subjected to abiotic stress conditions. Results are presented as 
mean relative expression with SD from three biological replicates after 
normalization using the best combination of reference genes recom- 
mended by geNorm and RefFinder (see Table 8) for abiotic stress 
conditions. Dehydration (mannitol), drought (no water), heat (42 C C) and 
cold (4"C) stresses are presented. Different letters on the bars indicate 
significant differences at the P<0.05 level as tested by Tukey's Range 
Test. 

doi:1 0.1 371 /journal.pone.01 06308.g004 

days post bombardment using qRT-PCR. Normalization with the 
recommended reference genes (PP2A, TIP41 and UBC2) showed 
similar effects on the relative expression patterns of gas, gfp and 
hpt transgenes in the calli of all three genotypes (Figure 5) as 
observed for PEPC in leaves (Figure 2), whereas the combination 
of the two {PP2A+TIP41) and the three (PP2A+TIP41+UBC2) 
reference genes exhibited more reliable transcript quantification. 
In general, expression analyses revealed that relative quantification 
of all three transgenes were higher in pearl millet genotype 
ICMT01004 and IPCI1466 compared to IP300088 (Figure 5). 

Discussion 

Transcriptome changes occurring during developmental pro- 
cesses and/or adverse environmental conditions are experiencing 
a growing research interest to understand the gene regulatory 
networks that control agronomically and economically important 
traits e.g. enhanced crop yield and biomass production under high 
atmospheric C0 2 or abiotic stress in the Panicoideae grasses 
including pearl millet. Transcriptomics data from microarray and 
next generation sequencing analyses should be validated using 
qRT-PCR [41]. QRT-PCR provides a useful tool to study 



OPP2A 

□ TIP41 

MUBC2 

MPP2A+TIP41 

MPP2A+TIP41+UBC2 




ICMR01004 



IPCI1466 



IP300088 



Figure 5. Validation of expression of gus, gfp and hpt 
transgenes using optimal number of control reference genes 
in hygromycin resistant calli from genotypes (A) ICMR01004, 
(B) IPCI1466 and (C) IP300088 after 30 days post particle 
bombardment-mediated transformation using pCAMBIA1201 
and pCAMBIAI 302, respectively. Results are presented as mean 
relative expression with SD from three biological replicates after 
normalization using the best combination of reference genes recom- 
mended by geNorm and RefFinder (see Table 8) for developmental 
tissue samples. Different letters on the bars indicate significant 
differences at the P<0.05 level as tested by Tukey's Range Test. 
doi:1 0.1 371/joumal.pone.01 06308.g005 

transcriptome changes in pearl millet because no genome 
sequence or microarray chip is available. Moreover, reliable 
transcript measurements using qRT-PCR analysis require accu- 
rate normalization against an appropriate internal control 
reference gene [9,28]. Normalization is important to adjust the 
variation introduced by various steps involved in the qRT-PCR 
such as quantity and quality of RNA samples, cDNAs, fluorescent 
fluctuations, sample-to-sample and/or well-to-well volume varia- 
tions [10]. Therefore, pearl millet requires an assessment of 
appropriate reference genes for accurate transcript normalization 
in gene expression studies using qRT-PCR. 

In this study, we demonstrated a comprehensive analysis of 18 
potential candidate reference genes which included both tradi- 
tional housekeeping genes like^CT, eEFla, GAPDH, TUA, UBC 
and UBQ5 and new candidate reference genes e.g. PEPKR, 
PP2A, TIP41 on 234 samples from developmental tissues, 
hormone treatments and abiotic stress conditions of three pearl 
millet genotypes. We carried out simple total RNA extraction 
protocols using the guanidinium thiocyanate-based kit [42] , which 



PLOS ONE | www.plosone.org 



13 



August 2014 | Volume 9 | Issue 8 | e106308 



Assessing Reference Genes for qRT-PCR in Pearl Millet 



yielded acceptable RNA quality and quantity from all samples 
including roots and seeds of three pearl millet genotypes as 
mentioned by the golden rules of qRT-PCR [11]. Previously 
published protocols using same guanidinium thiocyanate-based kit 
demonstrated satisfactory amount of high quality RNA from rice 
[43]. Since DNA contamination can result in inaccurate 
quantification of RNA abundance [44], we conducted a second 
gDNA wipeout reaction on the isolated RNAs after on-column 
DNase treatment following manufacturer recommended protocol 
to completely eliminate the detectable genomic DNA contamina- 
tion as verified by qRT-PCR for absence of any non-specific 
amplification. We primed cDNA synthesis using an optimized 
blend of oligo-dT and random primers to preferentially amplify 
the lowly abundant transcripts such as CYC in this study. In 
support of our finding, a weak expression of CYC was observed in 
rice [43]. However, the abundance of 25 S rRNA in different 
tissue, physiological conditions and pearl millet genotypes 
suggested the use of random hexamers to prime the reverse 
transcription reaction in our study. 

We performed the two-step qRT-PCR method to reduce the 
unwanted primer dimer formation using SYBR Green detection 
dye [8]. This method was also followed for the large scale 
expression profiling of transcription factors in rice [43] . Specific 
amplification with expected amplicon size of each primer pair 
from the RT-PCR was confirmed by agarose gel electrophoresis 
(Figure SI). In addition, the single peak melting curves in the qRT- 
PCR with no amplicon peak in the NTCs proved the absence of 
primer dimers or non-specific products (Figure S2). The PCR 
efficiency of each primer pair was calculated from the raw 
amplification curves (absolute fluorescence data) captured during 
the exponential phase of amplification of each qRT-PCR reaction 
using LinRegPCR [37]. Except for CYC and 25 S rRNA, which 
showed an average efficiency of 1.87±0.03 and 1.86±0.04, all the 
candidate reference genes exhibited mean efficiency values greater 
than 1.90 (Table SI) suggesting specific transcripts being amplified 
at least at 90% efficiency per cycle in the qRT-PCR reactions [45] . 
An identical range of PCR efficiencies were reported for many 
orthologous of selected candidate reference genes from Arabi- 
dopsis [12], rice [43] and common bean [46]. In this study the 
average Ct (Tables 1-2 and S3-S5) values of candidate reference 
genes varied within the recommended range of 22.8±3.1 to 
31.5±3.0 by qRT-PCR [47], except for 25S rRNA which showed 
Ct of 9. 1 ± 1.8 (Table 1). In support of our results, a low Ct 
(average Ct value of 8) of 25S rRNA gene was also observed in 
rice [13]. The ACt was calculated using the previously published 
method [39] and the precision of the assay was assessed using the 
CV. In general, our candidate reference genes showed CV<5% of 
Ct values, suggesting higher stability in expression levels under all 
experimental conditions. Therefore, our data demonstrated that 
the selected reference genes in this study are potential candidates 
for accurate normalization of gene expression by qRT-PCR after 
proper validation. In conjunction of our study, low CV<5% of Cq 
values of reference genes under abiotic stress conditions in 
common bean was also reported [46]. 

It has been suggested that the selection of optimal number of 
reference genes must be experimentally determined [48] . Howev- 
er, no single reference gene was found to have a stable expression 
under different experimental conditions [10,33] and nor a single 
method is enough to test for the stability of the candidate reference 
genes [31,33]. We used the algorithms executed by six different 
programs for proper stability ranking of the candidate reference 
genes. The SI [28] and ACt [29] methods calculate the variation 
of Ct and ACt values in pairwise genes, whereas BestKeeper 
estimates the variation in Ct values and reference genes showing 



SD<1 are considered the most stable [30]. However, the 
NormFinder [32] and geNorm [31] statistical algorithms allowed 
us to determine the stability ranking by calculating the SV and 
MV of each reference gene, respectively (Tables 3-7). In our study 
geNorm analyses revealed MV<1.5 for most of the genes under 
different experimental conditions (Tables 3-7), suggesting the 
potential stability of reference genes [31]. However, in the total 
experimental set PEPKR was the first ranked candidate gene by 
SI and BestKeeper, but ranked third by geNorm (Table 3); this 
could be due to the sensitivity of geNorm to the co-regulation of 
genes with similar expression patterns. In addition, geNorm is less 
affected by expression intensity of the reference genes [49] and 
allowed us to determine the optimal number of genes required to 
accurately normalize qRT-PCR data based on the V values [31]. 
We applied RefFinder [33] for recommended comprehensive 
ranking by combining all five above programs. Earlier reports on 
bamboo [17], strawberry [49] and leafy spurge [50] showed that 
these computational programs did not place the top ranked genes 
in identical order. According to our analysis, the six statistical 
programs ranked the candidate reference genes in various orders 
from best to worst, which could be due to different algorithm used 
by each program. Overall, new reference genes ranked better than 
the traditional housekeeping genes by most of the programs 
(Tables 3-8). Normalization using multiple reference genes is 
critical not only to obtain reliable gene expression results since 
normalization using single gene can be erroneous [9], but it also 
evaluates the expression stability of the selected reference genes 
during qRT-PCR. The geNorm analyses allowed us to identify 
optimal number of reference genes (Table 8) required for accurate 
normalization by calculating the V values at the suggested cut-off 
range of 0.15 [31]. 

In this study all the six computational methods suggested that 
PP2A, TIP41, UBC2, UBQ5 and ACT are the top 5 superior 
reference genes for accurate transcript normalization in pearl 
millet under different experimental conditions (Table 8). None of 
the traditional housekeeping genes qualified as the best reference 
gene for transcript normalization in total tissue across all the five 
experimental sets of pearl millet. Moreover, only UBQ5 and^4CT 
were found to be suitable for hormone treated, stress conditions 
and genotypes of pearl millet (Table 8), respectively. This is 
because expression stability of many housekeeping genes vary- 
considerably owning to their involvement in the cellular metab- 
olism and functions [26]. In accordance to our study, ACT was 
one of the best reference genes in foxtail millet [20] . In addition, 
ACT was shown be a good candidate reference gene for 
normalization of transcript data in rice [43] and strawberry 
[49]. Moreover, UBQwas found to be a suitable reference gene in 
mustard [21], poplar [28] and rice [13]. In the current study, 18S 
rRNA, 25 S rRNA and SAMDc were consistendy categorized as 
unsuitable, perhaps due to their inconsistency in gene expression 
by all the six programs (Tables 3-7), thereby rendering them 
inappropriate to use as reference gene. Similarly, poor stability of 
18S rRNA under abiotic stress conditions was reported in foxtail 
millet [20] . In conjunction with rice the high expression of 25S 
rRNA in this study makes it inappropriate for normalization of 
weakly expressed genes [13]. We observed significant variation of 
SAMDc expression pattern, which has been shown recently to be a 
poor reference gene in switchgrass [25]. The CYC, eEFltx and 
eIF4a were listed as variable genes in many studies [24,43], 
thereby limiting their use as reference genes in pearl millet as well. 
We found GAPDH as an inappropriate reference gene, which was 
also ranked unsuitable for bamboo [17], brachypodium [14] and 
rice [13]. In our study, another traditional housekeeping gene 
UBC2 ranked the third best reference genes after two novel 
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candidate reference genes, PP2A and TIP 41 for normalization in 
developmental tissue samples. The UBC encodes an ubiquitin- 
conjugating enzyme E2 involved in protein degradation through 
ubiquitination reactions and performed best among the three 
traditional housekeeping genes in leafy spurge [50]. However, in 
the current study two novel candidate genes, PP2A and TIP41 
resulted as superior reference genes compared to traditional 
housekeeping genes tested under different experimental condi- 
tions. This finding is in agreement with previous reports where 
PP2A and TIP41 combination was most suitable for abiotic stress 
conditions in caragana [27]. Recent reports demonstrated that 
PP2A and TIP41 were the most recommended stable reference 
genes for transcript normalization in tissue samples of numerous 
plant species [17,19,27]. 

The suitability of these reference genes to conduct transcrip- 
tomics studies was assessed by monitoring the expression profiles 
of three endogenous genes and transgenes in both untransformed 
and genetically transformed pearl millet tissues. The PEPC 
encodes a ubiquitous cytosolic enzyme in higher plants which 
catalyzes the irreversible carboxylation of phosphoenolpyruvate 
(PEP) to oxaloacetate (OAA), a four carbon compound, in the 
initial fixation of atmospheric C0 2 during C4 photosynthesis [51]. 
We noticed that transcript levels of PEPC were high in the flag 
leaf compared to nodal tissue in all the pearl millet genotypes 
studied (Figure 2). The ERF and DREB are AP2 binding 
transcription factors which regulate plant responses to several 
environmental stress conditions [52] and up-regulated under 
abiotic stresses [52] and hormone signaling [53], respectively. 
Transcript abundance of ERF illustrated differential expression 
pattern after accurate quantification using TIP41 and UBQ5 
under different hormone stimuli conditions (Figure 3). Currently, 
several reports have validated the optimum relative expression of 
DREB using appropriate reference genes under abiotic stress 
conditions [21,27]. In agreement with previous reports, we found 
DREB expression was up-regulated many fold in drought and heat 
stress conditions after accurate normalization using combination 
of reference genes (Figure 4). In addition, we provided evidence 
that these set of reference genes are also useful for transcript 
quantification in transformed pearl millet tissues, while incorpo- 
ration of multiple reference genes provides the most reliable 
expression pattern after precise normalization. 

Conclusions 

To the best of our knowledge this is the first comprehensive 
assessment of appropriate reference genes for accurate transcript 
normalization using qRT-PCR analyses in pearl millet. Stability 
ranking using computer based Stability Index, ACt, BestKeeper, 
NormFinder, geNorm and RefFinder programs recommended 
TIP41, PP2A, UBC2, UBQ5 and ACT as the best reference 
genes out of 18 potential candidate genes tested on different 
developmental and experimental conditions. This work will 
facilitate the developmental gene expression studies on C4 
photosynthesis and hormone cross-talk during abiotic stress 
conditions in pearl millet, a crop with limited genomic and 
transcriptomics information, and also benefit the scientific 
community for conducting experiments on related bioenergy crop 
species. 
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